Wavefield Propagator for Tilted Orthorhombic Media

ABSTRACT

Systems and methods that include receiving reservoir data of a hydrocarbon reservoir, receive an indication related to selection of a wavefield propagator, application of the wavefield propagator utilizing Fourier Finite Transforms and Finite Differences to model a wavefield associated with a Tilted Orthorhombic media representative of a region of a subsurface comprising the hydrocarbon reservoir, and processing the reservoir data in conjunction the wavefield propagator to generate an output for use with seismic exploration above a region of a subsurface comprising the hydrocarbon reservoir and containing structural or stratigraphic features conducive to a presence, migration, or accumulation of hydrocarbons.

BACKGROUND

The present disclosure relates generally to analyzing seismic data, andmore specifically, to utilizing a wavefield propagator to processseismic measurements.

This section is intended to introduce the reader to various aspects ofart that may be related to various aspects of the present disclosure,which are described and/or claimed below. This discussion is believed tobe helpful in providing the reader with background information tofacilitate a better understanding of the various aspects of the presentdisclosure. Accordingly, it should be understood that these statementsare to be read in this light, and not as admissions of prior art.

A seismic survey includes generating an image or map of a subsurfaceregion of the Earth by sending sound energy down into the ground andrecording the reflected sound energy that returns from the geologicallayers within the subsurface region. During a seismic survey, an energysource is placed at various locations on or above the surface region ofthe Earth, which may include hydrocarbon deposits. Each time the sourceis activated, the source generates a seismic (e.g., sound wave) signalthat travels downward through the Earth, is reflected, and, upon itsreturn, is recorded using one or more receivers disposed on or above thesubsurface region of the Earth. The seismic data recorded by thereceivers may then be used to create an image or profile of thecorresponding subsurface region.

SUMMARY

A summary of certain embodiments disclosed herein is set forth below. Itshould be understood that these aspects are presented merely to providethe reader with a brief summary of these certain embodiments and thatthese aspects are not intended to limit the scope of this disclosure.Indeed, this disclosure may encompass a variety of aspects that may notbe set forth below.

Seismic processing may include, for example, seismic data migration,full waveform inversion processing, reverse time migration, and othertechniques that utilize accurate (and often complex) representations ofa subsurface. The representation of the subsurface that is generated maybe termed a wavefield propagator (or a wave propagator) and varioustechniques may be applied as wave propagation methods to generatedifferent types of wavefield propagators for use in seismic processing.In some embodiments, these propagation methods may include apseudo-spectral (PS) wave propagation method, a tilted transverselyisotropic (TTI) wave propagation method (such as a TTI Fourier FiniteDifferences (FFD) wave propagation method), and a Tilted OrthorhombicFFD (TORFFD) wave propagation method (which, for example, may be a chainoperator of Fourier Finite Transforms and Finite Differences). As willbe discussed in greater detail below, use of the TORFFD wave propagationmethod may allow for greater resolution images of a region of asubsurface comprising a hydrocarbon reservoir and containing structuralor stratigraphic features conducive to a presence, migration, oraccumulation of hydrocarbons generated via processing of receivedseismic data relative to a TTI wave propagation method, while improvingat least on processing time required utilizing a pseudo-spectral wavepropagation method.

BRIEF DESCRIPTION OF THE DRAWINGS

Various aspects of this disclosure may be better understood upon readingthe following detailed description and upon reference to the drawings inwhich:

FIG. 1 is a flow chart of various processes that may be performed basedon analysis of seismic data acquired via a seismic survey system, inaccordance with embodiments presented herein;

FIG. 2 is a schematic diagram of a marine survey system in a marineenvironment, in accordance with embodiments presented herein;

FIG. 3 is a schematic diagram of a second marine survey system in amarine environment, in accordance with embodiments presented herein

FIG. 4 is a block diagram of a computing system that may performoperations described herein based on data acquired via the marine surveysystem of FIG. 2 and/or the second marine survey system of FIG. 3, inaccordance with embodiments presented herein;

FIG. 5 is a flow chart of a method for generating a seismic image viathe computing system of FIG. 4, in accordance with embodiments presentedherein;

FIG. 6 illustrates pseudo-acoustic wavefield snapshots for two TiltedOrthorhombic modeling techniques applied by the computing system of FIG.4, in accordance with embodiments presented herein;

FIG. 7 illustrates an example of a first seismic image generated by thecomputing system of FIG. 4 utilizing a TTI wave propagation method, inaccordance with embodiments presented herein;

FIG. 8 illustrates an example of a first seismic image generated by thecomputing system of FIG. 4 utilizing a pseudo-spectral wave propagationmethod on Tilted Orthorhombic media, in accordance with embodimentspresented herein; and

FIG. 9 illustrates an example of a first seismic image generated by thecomputing system of FIG. 4 utilizing a FFD wave propagation method onTilted Orthorhombic media, in accordance with embodiments presentedherein.

DETAILED DESCRIPTION OF SPECIFIC EMBODIMENTS

One or more specific embodiments will be described below. In an effortto provide a concise description of these embodiments, not all featuresof an actual implementation are described in the specification. Itshould be appreciated that in the development of any such actualimplementation, as in any engineering or design project, numerousimplementation-specific decisions must be made to achieve thedevelopers' specific goals, such as compliance with system-related andbusiness-related constraints, which may vary from one implementation toanother. Moreover, it should be appreciated that such a developmenteffort might be complex and time consuming, but would nevertheless be aroutine undertaking of design, fabrication, and manufacture for those ofordinary skill having the benefit of this disclosure.

Seismic data may provide valuable information with regard to thedescription such as the location and/or change of hydrocarbon depositswithin a subsurface region of the Earth. Different processingtechniques, such as the application of different wave propagationmethods for particularly selected media (e.g., for selectedrepresentations of a subsurface region of the Earth as wavefieldpropagators) result in different images generated by the processing ofthe seismic data and/or different amounts of time to generate thedifferent images generated by the processing of the seismic data.Accordingly, there may be situations in which a particular chosen media(or wavefield propagator) may benefit from a particular type of wavepropagation method. The present application details at least one wavepropagation method that may beneficially be used in conjunction withTilted Orthorhombic (TOR) media.

By way of introduction, seismic data may be acquired using a variety ofseismic survey systems and techniques, two of which are discussed withrespect to FIG. 2 and FIG. 3. Regardless of the seismic data gatheringtechnique utilized, after the seismic data is acquired, a computingsystem may analyze the acquired seismic data and may use the results ofthe seismic data analysis (e.g., seismogram, map of geologicalformations, etc.) to perform various operations within the hydrocarbonexploration and production industries. For instance, FIG. 1 illustratesa flow chart of a method 10 that details various processes that may beundertaken based on the analysis of the acquired seismic data. Althoughthe method 10 is described in a particular order, it should be notedthat the method 10 may be performed in any suitable order.

Referring now to FIG. 1, at block 12, locations and properties ofhydrocarbon deposits within a subsurface region of the Earth associatedwith the respective seismic survey may be determined based on theanalyzed seismic data. In one embodiment, the seismic data acquired maybe analyzed to generate a map or profile that illustrates variousgeological formations within the subsurface region. Based on theidentified locations and properties of the hydrocarbon deposits, atblock 14, certain positions or parts of the subsurface region may beexplored. That is, hydrocarbon exploration organizations may use thelocations of the hydrocarbon deposits to determine locations at thesurface of the subsurface region to drill into the Earth. As such, thehydrocarbon exploration organizations may use the locations andproperties of the hydrocarbon deposits and the associated overburdens todetermine a path along which to drill into the Earth, how to drill intothe Earth, and the like.

After exploration equipment has been placed within the subsurfaceregion, at block 16, the hydrocarbons that are stored in the hydrocarbondeposits may be produced via natural flowing wells, artificial liftwells, and the like. At block 18, the produced hydrocarbons may betransported to refineries and the like via transport vehicles,pipelines, and the like. At block 20, the produced hydrocarbons may beprocessed according to various refining procedures to develop differentproducts using the hydrocarbons.

It should be noted that the processes discussed with regard to themethod 10 may include other suitable processes that may be based on thelocations and properties of hydrocarbon deposits as indicated in theseismic data acquired via one or more seismic survey. As such, it shouldbe understood that the processes described above are not intended todepict an exhaustive list of processes that may be performed afterdetermining the locations and properties of hydrocarbon deposits withinthe subsurface region.

With the foregoing in mind, FIG. 2 is a schematic diagram of a marinesurvey system 22 (e.g., for use in conjunction with block 12 of FIG. 1)that may be employed to acquire seismic data (e.g., waveforms) regardinga subsurface region of the Earth in a marine environment. Generally, amarine seismic survey using the marine survey system 22 may be conductedin an ocean 24 or other body of water over a subsurface region 26 of theEarth that lies beneath a seafloor 28.

The marine survey system 22 may include a vessel 30, one or more seismicsources 32, a (seismic) streamer 34, one or more (seismic) receivers 36,and/or other equipment that may assist in acquiring seismic imagesrepresentative of geological formations within a subsurface region 26 ofthe Earth. The vessel 30 may tow the seismic source(s) 32 (e.g., an airgun array) that may produce energy, such as sound waves (e.g., seismicwaveforms), that is directed at a seafloor 28. The vessel 30 may alsotow the streamer 34 having a receiver 36 (e.g., hydrophones) that mayacquire seismic waveforms that represent the energy output by theseismic source(s) 32 subsequent to being reflected off of variousgeological formations (e.g., salt domes, faults, folds, etc.) within thesubsurface region 26. Additionally, although the description of themarine survey system 22 is described with one seismic source 32(represented in FIG. 2 as an air gun array) and one receiver 36(represented in FIG. 2 as a set of hydrophones), it should be noted thatthe marine survey system 22 may include multiple seismic sources 32 andmultiple receivers 36. In the same manner, although the abovedescriptions of the marine survey system 22 is described with onestreamer 34, it should be noted that the marine survey system 22 mayinclude multiple streamers similar to streamer 34. In addition,additional vessels 30 may include additional source(s) 32, streamer(s)34, and the like to perform the operations of the marine survey system22.

FIG. 3 is a block diagram of a land survey system 38 (e.g., for use inconjunction with block 12 of FIG. 1) that may be employed to obtaininformation regarding the subsurface region 26 of the Earth in anon-marine environment. The land survey system 38 may include aland-based seismic source 40 and land-based receiver 44. In someembodiments, the land survey system 38 may include multiple land-basedseismic sources 40 and one or more land-based receivers 44 and 46.Indeed, for discussion purposes, the land survey system 38 includes aland-based seismic source 40 and two land-based receivers 44 and 46. Theland-based seismic source 40 (e.g., seismic vibrator) that may bedisposed on a surface 42 of the Earth above the subsurface region 26 ofinterest. The land-based seismic source 40 may produce energy (e.g.,sound waves, seismic waveforms) that is directed at the subsurfaceregion 26 of the Earth. Upon reaching various geological formations(e.g., salt domes, faults, folds) within the subsurface region 26 theenergy output by the land-based seismic source 40 may be reflected offof the geological formations and acquired or recorded by one or moreland-based receivers (e.g., 44 and 46).

In some embodiments, the land-based receivers 44 and 46 may be dispersedacross the surface 42 of the Earth to form a grid-like pattern. As such,each land-based receiver 44 or 46 may receive a reflected seismicwaveform in response to energy being directed at the subsurface region26 via the seismic source 40. In some cases, one seismic waveformproduced by the seismic source 40 may be reflected off of differentgeological formations and received by different receivers. For example,as shown in FIG. 3, the seismic source 40 may output energy that may bedirected at the subsurface region 26 as seismic waveform 48. A firstreceiver 44 may receive the reflection of the seismic waveform 48 off ofone geological formation and a second receiver 46 may receive thereflection of the seismic waveform 48 off of a different geologicalformation. As such, the first receiver 44 may receive a reflectedseismic waveform 50 and the second receiver 46 may receive a reflectedseismic waveform 52.

Regardless of how the seismic data is acquired, a computing system(e.g., for use in conjunction with block 12 of FIG. 1) may analyze theseismic waveforms acquired by the receivers 36, 44, 46 to determineseismic information regarding the geological structure, the location andproperty of hydrocarbon deposits, and the like within the subsurfaceregion 26. FIG. 4 is a block diagram of an example of such a computingsystem 60 that may perform various data analysis operations to analyzethe seismic data acquired by the receivers 36, 44, 46 to determine thestructure and/or predict seismic properties of the geological formationswithin the subsurface region 26.

Referring now to FIG. 4, the computing system 60 may include acommunication component 62, a processor 64, memory 66, storage 68,input/output (I/O) ports 70, and a display 72. In some embodiments, thecomputing system 60 may omit one or more of the display 72, thecommunication component 62, and/or the input/output (I/O) ports 70. Thecommunication component 62 may be a wireless or wired communicationcomponent that may facilitate communication between the receivers 36,44, 46, one or more databases 74, other computing devices, and/or othercommunication capable devices. In one embodiment, the computing system60 may receive receiver data 76 (e.g., seismic data, seismograms, etc.)via a network component, the database 74, or the like. The processor 64of the computing system 60 may analyze or process the receiver data 76to ascertain various features regarding geological formations within thesubsurface region 26 of the Earth.

The processor 64 may be any type of computer processor or microprocessorcapable of executing computer-executable code. The processor 64 may alsoinclude multiple processors that may perform the operations describedbelow. The memory 66 and the storage 68 may be any suitable articles ofmanufacture that can serve as media to store processor-executable code,data, or the like. These articles of manufacture may representcomputer-readable media (e.g., any suitable form of memory or storage)that may store the processor-executable code used by the processor 64 toperform the presently disclosed techniques. Generally, the processor 64may execute software applications that include programs that processseismic data acquired via receivers of a seismic survey according to theembodiments described herein.

The memory 66 and the storage 68 may also be used to store the data,analysis of the data, the software applications, and the like. Thememory 66 and the storage 68 may represent non-transitorycomputer-readable media (e.g., any suitable form of memory or storage)that may store the processor-executable code used by the processor 64 toperform various techniques described herein. It should be noted thatnon-transitory merely indicates that the media is tangible and not asignal.

The I/O ports 70 may be interfaces that may couple to other peripheralcomponents such as input devices (e.g., keyboard, mouse), sensors,input/output (I/O) modules, and the like. I/O ports 70 may enable thecomputing system 60 to communicate with the other devices in the marinesurvey system 22, the land survey system 38, or the like via the I/Oports 70.

The display 72 may depict visualizations associated with software orexecutable code being processed by the processor 64. In one embodiment,the display 72 may be a touch display capable of receiving inputs from auser of the computing system 60. The display 72 may also be used to viewand analyze results of the analysis of the acquired seismic data todetermine the geological formations within the subsurface region 26, thelocation and property of hydrocarbon deposits within the subsurfaceregion 26, predictions of seismic properties associated with one or morewells in the subsurface region 26, and the like. The display 72 may beany suitable type of display, such as a liquid crystal display (LCD),plasma display, or an organic light emitting diode (OLED) display, forexample. In addition to depicting the visualization described herein viathe display 72, it should be noted that the computing system 60 may alsodepict the visualization via other tangible elements, such as paper(e.g., via printing) and the like.

With the foregoing in mind, the present techniques described herein mayalso be performed using a supercomputer that employs multiple computingsystems 60, a cloud-computing system, or the like to distributeprocesses to be performed across multiple computing systems 60. In thiscase, each computing system 60 operating as part of a super computer maynot include each component listed as part of the computing system 60.For example, each computing system 60 may not include the display 72since multiple displays 72 may not be useful to for a supercomputerdesigned to continuously process seismic data.

After performing various types of seismic data processing, the computingsystem 60 may store the results of the analysis in one or more databases74. The databases 74 may be communicatively coupled to a network thatmay transmit and receive data to and from the computing system 60 viathe communication component 62. In addition, the databases 74 may storeinformation regarding the subsurface region 26, such as previousseismograms, geological sample data, seismic images, and the likeregarding the subsurface region 26.

Although the components described above have been discussed with regardto the computing system 60, it should be noted that similar componentsmay make up the computing system 60. Moreover, the computing system 60may also be part of the marine survey system 22 or the land surveysystem 38, and thus may monitor and control certain operations of thesources 32 or 40, the receivers 36, 44, 46, and the like. Further, itshould be noted that the listed components are provided as examplecomponents and the embodiments described herein are not to be limited tothe components described with reference to FIG. 4.

In some embodiments, the computing system 60 may generate atwo-dimensional representation or a three-dimensional representation ofthe subsurface region 26 based on the seismic data received via thereceivers mentioned above. Additionally, seismic data associated withmultiple source/receiver combinations may be combined to create a nearcontinuous profile of the subsurface region 26 that can extend for somedistance. In a two-dimensional (2-D) seismic survey, the receiverlocations may be placed along a single line, whereas in athree-dimensional (3-D) survey the receiver locations may be distributedacross the surface in a grid pattern. As such, a 2-D seismic survey mayprovide a cross sectional picture (vertical slice) of the Earth layersas they exist directly beneath the recording locations. A 3-D seismicsurvey, on the other hand, may create a data “cube” or volume that maycorrespond to a 3-D picture of the subsurface region 26.

In addition, a 4-D (or time-lapse) seismic survey may include seismicdata acquired during a 3-D survey at multiple times. Using the differentseismic images acquired at different times, the computing system 60 maycompare the two images to identify changes in the subsurface region 26.

In any case, a seismic survey may be composed of a very large number ofindividual seismic recordings or traces. As such, the computing system60 may be employed to analyze the acquired seismic data to obtain animage representative of the subsurface region 26 and to determinelocations and properties of hydrocarbon deposits. To that end, a varietyof seismic data processing algorithms may be used to remove noise fromthe acquired seismic data, migrate the pre-processed seismic data,identify shifts between multiple seismic images, align multiple seismicimages, and the like.

After the computing system 60 analyzes the acquired seismic data, theresults of the seismic data analysis (e.g., seismogram, seismic images,map of geological formations, etc.) may be used to perform variousoperations within the hydrocarbon exploration and production industries.For instance, as described above, the acquired seismic data may be usedto perform the method 10 of FIG. 1 that details various processes thatmay be undertaken based on the analysis of the acquired seismic data.

In some embodiments, the results of the seismic data analysis may begenerated in conjunction with a seismic processing scheme such as, forexample, the method 78 illustrated in FIG. 5. As illustrated, method 78includes a seismic processing sequence that includes seismic datacollection in step 80, editing of the seismic data in step 82, initialprocessing in step 84, and signal processing, conditioning, and imaging(which may, for example, include production of imaged sections orvolumes as well as the wave propagation methods described herein used,for example, in connection with the illustrated migration process(es))in step 86 prior to any interpretation of the seismic data, any furtherimage enhancement consistent with the exploration objectives desired,generation of attributes from the processed seismic data,reinterpretation of the seismic data as needed, and determination and/orgeneration of a drilling prospect or other seismic survey applications.As a result of the processing method 78, location of hydrocarbons withina subsurface region 26 may be identified. Techniques for seismicmodeling generally, and more specifically, wave propagation methods thatmay beneficially be used in conjunction with Tilted Orthorhombic (TOR)media to improve upon the processing results of seismic data so as toallow for better detection of hydrocarbons within a subsurface region26.

In one embodiment, a two-way scalar wavefield propagator for TOR mediais generated. The wavefield propagator for TOR media may adopt apseudo-acoustic P-wave dispersion relation in the corresponding media.The wavefield propagator for TOR media combines Fast Fourier Transform(FFT) and Finite Differences (FD) operations to result in a TiltedOrthorhombic FFD (TORFFD) wave propagation method which, as a scalaroperator, can also overcome the coupling of qP-waves and qSV-wavespresent in the TOR media.

When used for TOR media, the TORFFD wave propagation method (andresultant wavefield propagator) may, for example, utilize two points perwavefield sampling in 3D space, similar to a pseudo-spectral (PS) wavepropagation method for TOR media. However, efficiencies are present byutilizing the TORFFD wave propagation method in place of a PS wavepropagation method, for example, utilizing the TORFFD wave propagationmethod includes the calculation of 2 forward FFTs, 3 inverse FFTs, and 1FD operation in each time step compared to 30 FFTs for the PS wavepropagation method. Furthermore, use of the TORFFD wave propagationmethod allows for a larger time step relative to a PS wave propagationmethod, for example, having a courant number being as large as 0.555 aswell as accelerated processing times (e.g. up to approximately two timesas fast, approximately three times as fast, approximately four times asfast, approximately five times as fast, or approximately six times asfast a PS wave propagation method) while achieving resultant accuracygains relative to a PS wave propagation method.

The theoretical underpinnings of the present TORFFD wave propagationmethod and associated wave propagator are described below. An acousticwave equation used in seismic modeling and reverse time migration (e.g.,a technique for imaging seismic data) is expressed as:

$\begin{matrix}{{\frac{\partial^{2}p}{\partial t^{2}} = {{v(x)}^{2}{\nabla^{2}p}}},} & \left( {{Equation}\mspace{20mu} 1} \right)\end{matrix}$

As set forth in Equation 1, p(x, t) is the seismic pressure wavefield,and v(x) is the propagation velocity. An acoustic pressure solution forEquation 1 may be represented as:

$\begin{matrix}{{{p\left( {x,{t + {\Delta \; t}}} \right)} + {p\left( {x,{t - {\Delta \; t}}} \right)} - {2{p\left( {x,t} \right)}}} = {2{\int_{- \infty}^{+ \infty}{{\hat{p}\left( {k,t} \right)}\left( {{\cos \left( {{k}v\; \Delta \; t} \right)} - 1} \right)e^{{- {ik}} \cdot x}{{dk}.}}}}} & \left( {{Equation}\mspace{20mu} 2} \right)\end{matrix}$

Equation 2 may operate as a solution in the case of a constant-velocitymedium with the aid of FFT. In a model with variable velocity, Equation2 can provide an approximation by replacing v with v(x). Further, in thecase of TOR, the term v(x) |k| on the right-hand side of Equation 2 canbe replaced with the acoustic approximation f (v, k), yielding:

$\begin{matrix}{{{f\left( {v,k} \right)} = \sqrt{{v_{zz}{\hat{k}}_{z}^{2}} + {v_{xx}{\hat{k}}_{x}^{2}} + {v_{yy}{\hat{k}}_{y}^{2}} - {v_{yz}\frac{{\hat{k}}_{y}^{2}{\hat{k}}_{z}^{2}}{{k}^{2}}} - {v_{xz}\frac{{\hat{k}}_{x}^{2}{\hat{k}}_{z}^{2}}{{k}^{2}}} - {v_{xy}\frac{{\hat{k}}_{x}^{2}{\hat{k}}_{y}^{2}}{{k}^{2}}}}},\begin{matrix}{{v = \left( {v_{z},ɛ_{1},ɛ_{2},\delta_{1},\delta_{2},\delta_{3},\theta,\varphi,\psi} \right)},} \\{{k = \left( {k_{x},k_{y},k_{z}} \right)},} \\{{v_{zz} = v_{z}^{2}},} \\{{v_{xx} = {v_{zz}\left( {1 + {2ɛ_{2}}} \right)}},} \\{{v_{yy} = {v_{zz}\left( {1 + {2ɛ_{1}}} \right)}},} \\{{v_{xy} = {2{v_{xx}\left( {\delta_{3} - \frac{ɛ_{1} - ɛ_{2}}{1 + {2ɛ_{2}}}} \right)}}},} \\{v_{xz} = {2{v_{zz}\left( {\delta_{2} + ɛ_{2}} \right)}}} \\{{v_{yz} = {2{v_{zz}\left( {\delta_{1} - ɛ_{1}} \right)}}},}\end{matrix}} & \left( {{Equation}\mspace{20mu} 3} \right)\end{matrix}$

As set forth in Equation 3, θ is the dip angle measured with respect tovertical, φ is the azimuth angle between the projection of the TIsymmetry axis in the horizontal plane and the original X-coordinate, ψis the rotation angle around the TI axis, ε₁, ε₂, δ₁, δ₂, and δ₃ areextended Thomsen parameters for orthorhombic media, v_(z) is the P-wavephase velocity along the TI symmetry axis, with the direction of (c_(x),c_(y), c_(z))={sin θ cos φ, sin θ sin φ, cos θ) while the “fast”velocity is along the direction of (r_(x), r_(y), r_(z))={cos θ cos φcos ψ−sin φ sin ψ, cos θ sin φ cos ψ+cos φ sin ψ, −sin θ cos ψ}, and{circumflex over (k)}_(x), {circumflex over (k)}_(y), and {circumflexover (k)}_(z) represent the wavenumbers evaluated in a rotatedcoordinate system aligned with a new axis.

When using Equation 2 to propagate waves in heterogeneous media, a FFTcan no longer be applied directly for the inverse Fourier transform fromthe wavenumber domain back to the space domain. As a solution, FourierFinite Differences (FFD), as a chain operator of the FFT and FDoperators may be utilized such that at each time step, the FFD methodfirst propagates the current wavefield using global referenceparameters, and then applies the wavefield correction according tovariations of local model parameters by Finite Differences. The presenttechniques for a wavefield propagator for TOR media may utilize asimilar approach by transforming the space-wavenumber mixed-domain termon the right-hand side of Equation 2 into a product of two terms:

$\begin{matrix}{{{{{\left\lbrack {{\cos \left( {{f\left( {v,k} \right)}\Delta \; t} \right)} - 1} \right\rbrack =}\quad}\left\lbrack \frac{{\cos \left( {{f\left( {v_{0},k} \right)}\Delta \; t} \right)} - 1}{{k}^{2}} \right\rbrack}\left\lbrack \frac{\left( {{\cos \left( {{f\left( {v,k} \right)}\Delta \; t} \right)} - 1} \right){k}^{2}}{{\cos \left( {{f\left( {v_{0},k} \right)}\Delta \; t} \right)} - 1} \right\rbrack},} & \left( {{Equation}\mspace{20mu} 4} \right)\end{matrix}$

As set forth in Equation 4, v₀ represents reference parameters. Thefirst term on the right hand side of Equation 4 is purely a function ofk. So we can use the inverse FFT to return to the space domain duringthe calculation. The remaining term G(v, k) may be solved bytransforming the phase shift in the wavenumber domain to finitedifferences in space with the following approximation:

$\begin{matrix}{{G\left( {v,k} \right)} = {\frac{\left( {{\cos \left( {{f\left( {v,k} \right)}\Delta \; t} \right)} - 1} \right){k}^{2}}{{\cos \left( {{f\left( {v_{0},k} \right)}\Delta \; t} \right)} - 1} \approx {a + {2{\sum\limits_{n = 1}^{3}\left\lbrack {{b_{n}{\cos \left( {k_{n}\Delta \; x_{n}} \right)}} + {d_{n}{\cos \left( {2k_{n}\Delta \; x_{n}} \right)}}} \right\rbrack}} + {2{\sum\limits_{n = 1}^{3}\left\lbrack {{{{{cp}_{n}{\cos \left( {{k_{i}\Delta \; x_{i}} + {k_{j}\Delta \; x_{j}}} \right)}} + \left. \quad{c\; m_{n}{\cos \left( {{k_{i}\Delta \; x_{i}} - {k_{j}\Delta \mspace{11mu} x_{j}}} \right)}} \right\rbrack} = {H\left( {v,k} \right)}},} \right.}}}}} & \left( {{Equation}\mspace{20mu} 5} \right)\end{matrix}$

As set forth in Equation 5, i, j=1, 2, 3; i≠j; i, j≠n and a, b_(n),c_(pn), c_(mn) and d_(n) are determined from a Taylor expansion aroundk=0, illustrated in Table 1 as variables of local model parameters inspace and serve as the 4^(th) order FD coefficients (with i=1, 2, 3):

TABLE 1 a${- 2}{\sum\limits_{n = 1}^{3}\left( {b_{n} + {cm}_{n} + {cp}_{n} + d_{n}} \right)}$b₁${{- 4}d_{1}} - {bb}_{2} - {bb}_{3} - \frac{g_{1}}{\Delta \; x^{2}}$b₂${{- 4}d_{2}} - {bb}_{3} - {bb}_{1} - \frac{g_{2}}{\Delta \; y^{2}}$b₃${{- 4}d_{3}} - {bb}_{1} - {bb}_{2} - \frac{g_{3}}{\Delta \; z^{2}}$d₁$\frac{1}{12}\left( {\frac{g_{1}}{\Delta \; x^{2}} + \frac{q_{1}}{\Delta \; x^{4}}} \right)$d₂$\frac{1}{12}\left( {\frac{g_{2}}{\Delta \; y^{2}} + \frac{q_{2}}{\Delta \; y^{4}}} \right)$d₃$\frac{1}{12}\left( {\frac{g_{3}}{\Delta \; z^{2}} + \frac{q_{3}}{\Delta \; z^{4}}} \right)$cm_(i) $\frac{1}{2}\left( {{bb}_{i} - {aa}_{i}} \right)$ cp_(i)$\frac{1}{2}\left( {{bb}_{i} + {aa}_{i}} \right)$

Table 2 provides further expressions of variables that appear in Table1, in which w_(0i) has a similar expression as w_(i) with all the spacevariables v replaced by the corresponding constant model parameters v₀.

TABLE 2 cs₁ c_(x)² cs₂ c_(y)² cs₃ c_(z)² cs₄$\frac{1}{2}\left( {c_{y} + c_{z}} \right)^{2}$ cs₅$\frac{1}{2}\left( {c_{x} + c_{z}} \right)^{2}$ cs₆$\frac{1}{2}\left( {c_{x} + c_{y}} \right)^{2}$ cs₇$\frac{1}{2}\left( {c_{y} - c_{z}} \right)^{2}$ cs₈$\frac{1}{2}\left( {c_{x} - c_{z}} \right)^{2}$ cs₉$\frac{1}{2}\left( {c_{x} - c_{y}} \right)^{2}$ rs₁ r_(x)² rs₂ r_(y)²rs₃ r_(z)² rs₄ $\frac{1}{2}\left( {r_{y} + r_{z}} \right)^{2}$ rs₅$\frac{1}{2}\left( {r_{x} + r_{z}} \right)^{2}$ rs₆$\frac{1}{2}\left( {r_{x} + r_{y}} \right)^{2}$ rs₇$\frac{1}{2}\left( {r_{y} - r_{z}} \right)^{2}$ rs₈$\frac{1}{2}\left( {r_{x} - r_{z}} \right)^{2}$ rs₉$\frac{1}{2}\left( {r_{x} - r_{y}} \right)^{2}$ w_(i) v_(zz)cs_(i) +v_(xx)rs_(i) + v_(yy)(1 − cs_(i) − rs_(i)) + v_(xz)cs_(i)rs_(i) +v_(yz)cs_(i)(1 − cs_(i) − rs_(i)) + v_(xy)rs_(i)(1 − cs_(i) − rs_(i))g_(i) $\frac{w_{i}}{w_{0i}}$ q_(i) Δt²w_(i)(1 − g_(i)) aa₁$\frac{g_{7} - g_{4}}{2\Delta \; y\; \Delta \; z}$ aa₂$\frac{g_{8} - g_{5}}{2\Delta \; z\; \Delta \; x}$ aa₃$\frac{g_{9} - g_{6}}{2\Delta \; x\; \Delta \; y}$ bb₁${\frac{{\Delta \; y^{2}} + {\Delta \; z^{2}}}{3\Delta \; y^{2}\Delta \; z^{2}}\left( {g_{4} + g_{7} - g_{2} - g_{3}} \right)} + \frac{{2q_{4}} + {2q_{7}} - q_{2} - q_{3}}{6\Delta \; y^{2}\Delta \; z^{2}}$bb₂${\frac{{\Delta \; z^{2}} + {\Delta \; x^{2}}}{3\Delta \; z^{2}\Delta \; x^{2}}\left( {g_{5} + g_{8} - g_{3} - g_{1}} \right)} + \frac{{2q_{5}} + {2q_{8}} - q_{3} - q_{1}}{6\Delta \; z^{2}\Delta \; x^{2}}$bb₃${\frac{{\Delta \; x^{2}} + {\Delta \; y^{2}}}{3\Delta \; x^{2}\Delta \; y^{2}}\left( {g_{6} + g_{9} - g_{1} - g_{2}} \right)} + \frac{{2q_{6}} + {2q_{9}} - q_{1} - q_{2}}{6\Delta \; x^{2}\Delta \; y^{2}}$

Extended to a case with TOR media, in order to get p(x, t+Δt) from p(x,t−Δt) and p(x, t), the following steps may be undertaken. Transformationof p(x, t) to {circumflex over (p)}(k, t) via a FFT may be undertaken.Next, {circumflex over (p)}(k, t) may be multiplied by

$\left\lbrack \frac{2\left( {{\cos \left( {{f\left( {v_{0},k} \right)}\Delta \; t} \right)} - 1} \right.}{k^{2}} \right\rbrack$

to get {circumflex over (q)}(k, t), which may then be transformed toq(k, t) through the use of inverse FFT. Finite differences may beapplied to q(k, t), utilizing the coefficients previously noted in Table1, to generate u(x, t). Transformation of u(x, t) to û(x, t) via a FFTmay be accomplished and correction terms u₀ in the wavenumber domain maybe applied, such that

${{\hat{u}}_{0}\left( {k,t} \right)} = {{{{\hat{u}}_{1}\left( {k,t} \right)}\frac{G\left( {v_{0},k} \right)}{H\left( {v_{0},k} \right)}\mspace{14mu} {and}\mspace{14mu} {{\hat{u}}_{1}\left( {k,t} \right)}} = {{{\hat{u}}_{1}\left( {k,t} \right)}\frac{G\left( {v_{1},k} \right)}{H\left( {v_{1},k} \right)}}}$

and, subsequently, transformation of û₀(k, t) and û₁(k, t) to u₀(k, t)and u₁(k, t), respectively, may be performed via the use of an inverseFFT. Subsequently, u₀(x, t) and u₁(x, t) may be combined, for example,according to local velocity variations to u(x, t) to result in p(x,t+Δt)←u(x, t+Δt)+2p(x, t)−p(x, t−Δt), whereby q(x, t), {circumflex over(q)}(k, t), u(x, t), û(k, t), u₀(x, t), û₀(k, t), u₁(x, t), and û₁(k,t), are temporary functions. Additionally, v₁ represents modelparameters at the location with minimum velocity, G(v₁, k) is the exactresponse for Equation 5, H(v₁, k) is an approximation by FD.Furthermore,

$\frac{G\left( {v_{1},k} \right)}{H\left( {v_{1},k} \right)}$

is used as a correction term for the wavefield in the wavenumber domainat locations with relatively low density while

$\frac{G\left( {v_{0},k} \right)}{H\left( {v_{0},k} \right)}$

is utilized as a second (e.g., the other) correction term.

The aforementioned process may be undertaken by the computing system 60as a TORFFD wave propagation method. In some embodiments, at least someof the aforementioned steps may be performed in an alternative order oromitted entirely. Moreover, while the steps may be appreciated as beingperformed by the processor 64 of the computing system 60, it should beunderstood that any suitable devices or systems, or combination ofsuitable devices or systems, including the processor 64, may perform theaforementioned steps, such as processing units or circuitry of computingdevices or systems external to, but communicatively coupled to, thecomputing system 60 and the execution of the steps may involve theprocessor 64 operating in conjunction with or utilizing software storedon a tangible machine readable medium to perform the aforementionedsteps in generating a TORFFD wave propagation method (and associatedwavefield propagator).

Additionally, it should be appreciated that the above described TORFFDTORFFD wave propagation method executed by the processor 64 represents aspecific improvement over conventional systems, and an improvedcomputing system 60 having improved functionality. In particular,utilizing the TORFFD wave propagation method and/or having a TORFFDwavefield propagator as part of the computing system improves, forexample, accuracy and efficiency with respect to alternate wavepropagation models and their associated waveform propagators by thecomputing system 60 when compared to performing seismic data analysisby, for example, computer systems performing conventional computeralgorithms and/or performing seismic data analysis without the TORFFDwave propagation method and/or its associated TORFFD wavefieldpropagator.

As a visual aid to represent some of the aforementioned advantages, FIG.6 illustrates examples of wave propagator responses in a TOR model. Apseudo-acoustic wavefield snapshot 88 in a constant TOR model by a PSwave propagation method is illustrated in FIG. 6. Likewise, a FFD method(i.e., the TORFFD wave propagation method) pseudo-acoustic wavefieldsnapshot 90 in a constant TOR model is additionally illustrated. As maybe observed, the pseudo-acoustic wavefield snapshot 90 possesses reduceddispersion than that present in pseudo-acoustic wavefield snapshot 88.Additionally, the process for generation of the pseudo-acousticwavefield snapshot 90 is up to six times faster than the process forgeneration of the pseudo-acoustic wavefield snapshot 88. Additionally,use of both the PS wave propagation method and the PS wave propagationmethod possess an advantage for pseudo-acoustic wave propagation in thatno coupling of qP-waves and qSv-waves is present.

As a second set of examples of the advantages of the TORFFD wavepropagation method, FIG. 7 illustrates a reverse time migration (RTM)image 92 generated utilizing a TTIFDD wave propagation method, FIG. 8illustrates an RTM image 94 generated utilizing a PS wave propagationmethod, and FIG. 5 illustrates an RTM image 96 generated utilizing theabove described the TORFFD wave propagation method. The runtime togenerate the RTM image 92 and the RTM image 96 was approximately 20minutes, while the runtime to generate the RTM image 94 wasapproximately two hours. Additionally, as may be appreciated, one canobserve some main events in RTM image 92; however, the fault planes inthe RTM image 92 do not appear to be clearly defined. In contrast, thefault boundary in the RTM image 96 is clearer than in RTM image 92,which allows for an improved ability to present strong azimuthalanisotropy. Additionally, events below or adjacent faults in the RTMimage 96 are more clearly imaged with respect to the RTM image 92. TheRTM image 94 likewise appears to provide improved imaging with respectto RTM image 92; however, it is noted that the runtime for thegeneration of RTM image 94 is approximately 6 times the runtime for thegeneration of the RTM image 96. Accordingly, utilization of the TORFFDwave propagation method appears to allow for the increased processingtime associated with a TTI FFD wave propagation method while (at least)maintaining or approximating the image quality in a generated RTM imagewith respect to a PS wave propagation method being utilized.

The specific embodiments described above have been shown by way ofexample, and it should be understood that these embodiments may besusceptible to various modifications and alternative forms. It should befurther understood that the claims are not intended to be limited to theparticular forms disclosed, but rather to cover all modifications,equivalents, and alternatives falling within the spirit and scope ofthis disclosure.

The techniques presented and claimed herein are referenced and appliedto material objects and concrete examples of a practical nature thatdemonstrably improve the present technical field and, as such, are notabstract, intangible or purely theoretical. Further, if any claimsappended to the end of this specification contain one or more elementsdesignated as “means for [perform]ing [a function] . . . ” or “step for[perform]ing [a function] . . . ”, it is intended that such elements areto be interpreted under 35 U.S.C. 112(f). However, for any claimscontaining elements designated in any other manner, it is intended thatsuch elements are not to be interpreted under 35 U.S.C. 112(f).

What is claimed is:
 1. A system, comprising: a processor, wherein theprocessor is configured to: receive reservoir data of a hydrocarbonreservoir; receive an indication related to selection of a wavefieldpropagator; apply the wavefield propagator utilizing Fourier FiniteTransforms and finite differences to model a wavefield associated with aTilted Orthorhombic media representative of a region of a subsurfacecomprising the hydrocarbon reservoir; and process the reservoir data inconjunction the wavefield propagator to generate an output for use withseismic exploration above a region of a subsurface comprising thehydrocarbon reservoir and containing structural or stratigraphicfeatures conducive to a presence, migration, or accumulation ofhydrocarbons.
 2. The system of claim 1, wherein the processor isconfigured to model the wavefield associated with the TiltedOrthorhombic media representative of the region of the subsurfacecomprising the hydrocarbon reservoir at least in part by performing aninitial transformation of a pressure attribute associated with thewavefield propagator into a transformed pressure attribute.
 3. Thesystem of claim 2, wherein the processor is configured to apply a FastFourier Transform to the pressure attribute to perform the initialtransformation of the pressure attribute associated with the wavefieldpropagator.
 4. The system of claim 2, wherein the processor isconfigured to model the wavefield associated with the TiltedOrthorhombic media representative of the region of the subsurfacecomprising the hydrocarbon reservoir at least in part by producing aproduct of the transformed pressure attribute with a propagationvelocity value to generate a first value.
 5. The system of claim 4,wherein the processor is configured to model the wavefield associatedwith the Tilted Orthorhombic media representative of the region of thesubsurface comprising the hydrocarbon reservoir at least in part by atleast in part by performing a transformation of the first value into atransformed first value.
 6. The system of claim 5, wherein the processoris configured to apply an inverse Fast Fourier Transform to the firstvalue to perform the transformation of the first value.
 7. The system ofclaim 5, wherein the processor is configured to model the wavefieldassociated with the Tilted Orthorhombic media representative of theregion of the subsurface comprising the hydrocarbon reservoir at leastin part by applying a finite differences operation utilizing apredetermined coefficient to the transformed first value to generate asecond value.
 8. The system of claim 7, wherein the processor isconfigured to model the wavefield associated with the TiltedOrthorhombic media representative of the region of the subsurfacecomprising the hydrocarbon reservoir at least in part by performing atransformation of the second value into a transformed second value. 9.The system of claim 8, wherein the processor is configured to apply aFast Fourier Transform to the second value to perform the transformationof the second value.
 10. The system of claim 8, wherein the processor isconfigured to model the wavefield associated with the TiltedOrthorhombic media representative of the region of the subsurfacecomprising the hydrocarbon reservoir at least in part by applying acorrection term in the wavenumber domain to the transformed second valueto generate a modified second value.
 11. The system of claim 10, whereinthe processor is configured to model the wavefield associated with theTilted Orthorhombic media representative of the region of the subsurfacecomprising the hydrocarbon reservoir at least in part by performing atransformation of the modified second value into a third value.
 12. Thesystem of claim 11, wherein the processor is configured to model thewavefield associated with the Tilted Orthorhombic media representativeof the region of the subsurface comprising the hydrocarbon reservoir atleast in part by applying a correction term in the wavenumber domain tothe second value to generate a second modified second value.
 13. Thesystem of claim 12, wherein the processor is configured to model thewavefield associated with the Tilted Orthorhombic media representativeof the region of the subsurface comprising the hydrocarbon reservoir atleast in part by performing a transformation of the second modifiedsecond value into a fourth value.
 14. The system of claim 13, whereinthe processor is configured to model the wavefield associated with theTilted Orthorhombic media representative of the region of the subsurfacecomprising the hydrocarbon reservoir at least in part by performing atransformation of the third value and the fourth value into atransformed third value and a transformed fourth value, respectively.15. The system of claim 14, wherein the processor is configured to applya Fast Fourier Transform to the third value and the fourth value toperform the transformation of the third value and the fourth value. 16.The system of claim 14, wherein the processor is configured to model thewavefield associated with the Tilted Orthorhombic media representativeof the region of the subsurface comprising the hydrocarbon reservoir atleast in part by combining the transformed third value with thetransformed fourth value.
 17. A method, comprising: applying a wavefieldpropagator utilizing Fourier Finite Transforms and Finite Differences tomodel a wavefield associated with a Tilted Orthorhombic mediarepresentative of a region of a subsurface comprising the hydrocarbonreservoir; and processing the reservoir data in conjunction thewavefield propagator to generate an output for use with seismicexploration above a region of a subsurface comprising a hydrocarbonreservoir and containing structural or stratigraphic features conduciveto a presence, migration, or accumulation of hydrocarbons.
 18. Themethod of claim 17, comprising modeling the wavefield associated withthe Tilted Orthorhombic media representative of the region of thesubsurface comprising the hydrocarbon reservoir at least in part by:performing an initial transformation of a pressure attribute associatedwith the wavefield propagator into a transformed pressure attribute;producing a product of the transformed pressure attribute with apropagation velocity value to generate a first value; performing atransformation of the first value into a transformed first value;applying a finite differences operation utilizing a predeterminedcoefficient to the transformed first value to generate a second value;performing a transformation of the second value into a transformedsecond value applying a correction term in the wavenumber domain to thetransformed second value to generate a modified second value; performinga transformation of the modified second value into a third value;applying a correction term in the wavenumber domain to the second valueto generate a second modified second value; performing a transformationof the second modified second value into a fourth value; performing atransformation of the third value and the fourth value into atransformed third value and a transformed fourth value, respectively;and combining the transformed third value with the transformed fourthvalue.
 19. One or more tangible, non-transitory, machine-readable mediacomprising instructions configured to cause a processor to: apply awavefield propagator utilizing Fourier Finite Transforms and FiniteDifferences to model a wavefield associated with a Tilted Orthorhombicmedia representative of a region of a subsurface comprising thehydrocarbon reservoir; and process the reservoir data in conjunction thewavefield propagator to generate an output for use with seismicexploration above a region of a subsurface comprising a hydrocarbonreservoir and containing structural or stratigraphic features conduciveto a presence, migration, or accumulation of hydrocarbons.
 20. The oneor more machine-readable media of claim 19, comprising instructionsconfigured to cause a processor to model the wavefield associated withthe Tilted Orthorhombic media representative of the region of thesubsurface comprising the hydrocarbon reservoir at least in part by:performing an initial transformation of a pressure attribute associatedwith the wavefield propagator into a transformed pressure attribute;producing a product of the transformed pressure attribute with apropagation velocity value to generate a first value; performing atransformation of the first value into a transformed first value;applying a finite differences operation utilizing a predeterminedcoefficient to the transformed first value to generate a second value;performing a transformation of the second value into a transformedsecond value applying a correction term in the wavenumber domain to thetransformed second value to generate a modified second value; performinga transformation of the modified second value into a third value;applying a correction term in the wavenumber domain to the second valueto generate a second modified second value; performing a transformationof the second modified second value into a fourth value; performing atransformation of the third value and the fourth value into atransformed third value and a transformed fourth value, respectively;and combining the transformed third value with the transformed fourthvalue.